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Prediction  of  Falling  Cylinder 
Through  Air-Water-Sediment 
Columns 

A  falling  rigid  body  through  air,  water,  and  sediment  is  investigated  experimentally  and 
theoretically.  Two  experiments  were  conducted  to  drop  rigid  cylinders  with  density  ratio 
around  1.8  into  shallow  water  (around  13  m  deep )  in  the  Monterey  Bay  (Exp-1)  and  into 
the  Naval  Postgraduate  School’s  swimming  pool  (Exp-2).  During  the  experiments,  we 
carefidly  observe  cylinder  track  and  burial  depth  while  simultaneously  taking  gravity 
cores  (in  Exp-1 ).  After  analyzing  the  gravity  cores,  we  obtain  the  bottom  sediment  density 
and  shear  strength  profiles.  The  theoretical  work  includes  the  development  of  a  3D  rigid 
body  impact  burial  prediction  model  (1MPACT35)  that  contains  three  components:  triple 
coordinate  transform  and  hydrodynamics  of  a  falling  rigid  object  in  a  single  medium  (air, 
water,  or  sediment )  and  in  multiple  media  (air-water  and  water-sediment  interfaces).  The 
model  predicts  the  rigid  body’s  trajectory  in  the  water  column  and  burial  depth  and 
orientation  in  the  sediment.  The  experimental  data  (burial  depth,  sediment  density,  and 
shear  strength)  show  the  capability  of  IMPACT3 5  in  predicting  the  cylinder’s  trajectory 
and  orientation  in  a  water  column  and  burial  depth  and  orientation  in  sediment. 
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1  Introduction 

Study  on  the  movement  of  a  rigid  body  in  fluid  has  wide  sci¬ 
entific  significance  and  technical  application.  The  scientific  stud¬ 
ies  of  the  hydrodynamics  of  a  rigid  cylinder  in  fluid  involve  the 
nonlinear  dynamics,  flight  theory,  body-fluid  interaction,  and  in¬ 
stability  theory.  The  body  forces  include  the  gravity  and  the  buoy¬ 
ancy  force.  The  hydrodynamic  forces  include  the  drag  and  lift 
forces  that  depend  on  the  fluid-to-body  velocity  and  the  impact 
force  as  the  body  penetrates  the  air-water  or  water-sediment  inter¬ 
faces.  Usually,  a  nonlinear  dynamical  system  is  needed  to  predict 
a  falling  rigid  body  in  fluid,  e.g.,  [1]. 

Recently,  the  scientific  problem  about  rigid  body  movement  in 
the  air-water-sediment  columns  drew  attention  to  the  naval  re¬ 
search.  This  is  due  to  the  threat  of  mines  in  the  naval  operations. 
Within  the  past  15  years  three  U.S.  ships,  the  USS  Samuel  B. 
Roberts  (FFG-58),  Tripoli  (LPF1-10),  and  Princeton  (CG-59)  have 
fallen  victim  to  mines.  Total  ship  damage  was  $125  million  while 
the  mines  cost  approximately  $30,000  [2].  Mines  have  evolved 
over  the  years  from  the  dumb  “horned”  contact  mines  that  dam¬ 
aged  the  Tripoli  and  Roberts  to  ones  that  are  relatively 
sophisticated — nonmagnetic  materials,  irregular  shapes,  anechoic 
coatings,  multiple  sensors,  and  ship  count  routines.  Despite  their 
increased  sophistication,  mines  remain  inexpensive  and  are  rela¬ 
tively  easy  to  manufacture,  keep,  and  place.  Water  mines  are  char¬ 
acterized  by  three  factors  [3,4]:  position  in  water  (bottom, 
moored,  rising,  and  floating),  method  of  delivery  (aircraft,  surface, 
and  subsurface),  and  method  of  actuation  (acoustic  and/or  mag¬ 
netic  influence,  pressure,  contact,  and  controlled).  Accurate  mine 
burial  predictions  are  inherently  difficult  to  make  because  of  un¬ 
certainties  in  both  mine  deployment  conditions  and  the  relevant 
environmental  parameters  [5],  The  U.S.  Navy  developed  opera- 
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tional  models  to  predict  the  environmental  parameters  for  mine 
burial  prediction  [6],  Recently,  statistical  methods  such  as  the 
Monte  Carlo  method  [7]  and  the  expert  system  method  [5]  have 
been  developed.  These  methods  need  a  core-physical  model  for 
describing  the  movement  of  falling  rigid  body  through  air-water- 
sediment  columns. 

When  the  rigid  body  is  cylindrical,  this  dynamical  system  can 
be  simplified  using  three  coordinate  systems:  earth-fixed  coordi¬ 
nate  (E-coordinate),  cylinder’s  main-axis  following  coordinate 
(M-coordinate),  and  hydrodynamic  force  following  coordinate  (F- 
coordinate).  The  origin  of  both  M-  and  F-coordinates  is  at  the 
cylinder’s  center  of  mass  (COM).  The  body  forces  and  their  mo¬ 
ments  are  easily  calculated  using  the  E-coordinate  system.  The 
hydrodynamic  forces  and  their  moments  are  easily  computed  us¬ 
ing  the  F-coordinate.  The  cylinder’s  moments  of  gyration  are  sim¬ 
ply  represented  using  the  M-coordinate.  Recently,  Chu  et  al.  [8] 
developed  a  recursive  model  to  predict  the  cylinder’s  translation 
velocity  and  orientation  in  the  water  column  (single  phase)  on  the 
base  of  the  triple  coordinate  transformation. 

To  extend  the  recursive  model  from  single  medium  (water  col¬ 
umn)  to  multi-media  (air,  water,  sediment),  a  falling  cylinder 
through  air-water  and  water-sediment  interfaces  (i.e.,  cylinder 
contacting  with  two  media)  should  be  particularly  analyzed.  The 
cylinder  is  decomposed  into  two  parts  with  each  one  contacting 
one  medium.  For  the  air-water  penetration,  the  cylinder  is  decom¬ 
posed  into  air  and  water  parts.  For  the  water-sediment  penetration, 
the  cylinder  is  decomposed  into  water  and  sediment  parts.  The 
body  forces  (such  as  the  buoyancy  force)  and  surface  forces  (such 
as  pressure  and  hydrodynamic  force)  are  computed  separately  for 
the  two  parts.  A  fully  three-dimensional  model  is  developed  for 
prediction  of  the  translation  velocity  and  orientation  of  falling 
rigid  cylinder  through  air,  water,  and  sediment.  Theoretical  model 
development  and  a  cylinder  drop  experiment  for  the  model  evalu¬ 
ation  are  depicted  in  this  paper. 

The  outline  of  this  paper  is  as  follows:  Section  2  depicts  the 
triple  coordinate  systems.  Section  3  describes  the  dynamics  for 
determining  the  cylinder’s  translation  velocity  and  orientation. 
Section  4  presents  the  equivalent  cylinder  method  for  computing 
hydrodynamic  forces  and  torques  when  the  cylinder  penetrates  the 
air-water  and  water-sediment  interfaces.  Section  5  describes 
forces  and  torques  in  air  and  water.  Section  6  describes  the  resis- 
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Fig.  1  M-coordinate  with  the  COM  as  the  origin  X  and  ( im,jm ) 
as  the  two  axes.  Here,  \  is  the  distance  between  the  COV  (B) 
and  COM  (X),  (/.,/?)  are  the  cylinder’s  length  and  radius. 


tance  from  sediments.  Section  7  shows  the  model  integration.  Sec¬ 
tion  8  describes  two  cylinder  drop  experiments  and  observational 
data  processing.  Section  9  presents  the  model-data  inter  compari¬ 
son.  The  conclusions  are  listed  in  Sec.  10. 


2  Triple  Coordinate  Systems 

Consider  an  axially  symmetric  cylinder  with  the  centers  of 
mass  (COM)  X  [or  called  gravity  center  (GC)  in  literatures]  and 
center  of  volume  (COV)  B  on  the  main  axis  (Fig.  1).  Let  ( L,R,x ) 
represent  the  cylinder’s  length,  radius,  and  the  distance  between 
the  two  points  (X,B).  The  positive  \  values  refer  to  nose-down 
case,  i.e.,  the  point  X  is  lower  than  the  point  B.  Three  coordinate 
systems  are  used  to  model  the  falling  cylinder  through  the  air, 
water,  and  sediment  phases:  earth-fixed  coordinate  (E-coordinate), 
main-axis  following  coordinate  (M-coordinate),  and  force  follow¬ 
ing  coordinate  (F-coordinate)  systems.  All  the  systems  are  three- 
dimensional,  orthogonal,  and  right-handed  [8]. 

2.1  E-Coordinate.  The  E-coordinate  is  represented  by 
FE(0,i,j,k)  with  the  origin  “O”  and  three  axes:  x,  y  axes  (hori¬ 
zontal)  with  the  unit  vectors  (i,j)  and  z  axis  (vertical)  with  the 
unit  vector  k  (upward  positive).  The  position  of  the  cylinder  is 
represented  by  the  position  of  the  COM, 


Fig.  2  Three  coordinate  systems 


X  =  xi  +  yj+zk,  (1) 

which  is  translation  of  the  cylinder.  The  translation  velocity  is 
given  by 

^/X 

—  =  V,  V  =  (u,v,w).  (2) 

at 

2.2  M-Coordinate.  Let  orientation  of  the  cylinder’s  main  axis 
(pointing  downward)  be  given  by  iM.  The  angle  between  iM  and  k 
is  denoted  by  i//2+tt/2.  Projection  of  the  vector  iM  onto  the  (x,y) 
plane  creates  angle  (if/ 3)  between  the  projection  and  the  x  axis 
(Fig.  2).  The  M-coordinate  is  represented  by  FM(X,iM,jM,kM) 
with  the  origin  “X”,  unit  vectors  (iM,jM,kM),  and  coordinates 
(*m.Xm,Zm)-  In  the  plane  consisting  of  vectors  iM  and  k  (passing 
through  the  point  M),  two  new  unit  vectors  (jM,kM)  are  defined 
with  jM  perpendicular  to  the  (iM-k)  plane,  and  kM  perpendicular 
to  iM  in  the  (iM,k)  plane.  The  unit  vectors  of  the  M-coordinate 
system  are  given  by  (Fig.  2) 

=  kM  =  iMXjM.  (3) 

The  M-coordinate  system  is  solely  determined  by  orientation  of 
the  cylinder’s  main  axis  iM.  Let  the  vector  P  be  represented  by  £P 
in  the  E-coordinate  and  by  A/P  in  the  M-coordinate,  and  let  be 
the  rotation  matrix  from  the  M-coordinate  to  the  E-coordinate, 
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Transformation  of  A/P  into  £P  contains  rotation  and  translation, 


4l’ 


^R(^3)MP  +  X.  (6) 

Let  the  cylinder  rotate  around  (iM,jM,kM)  with  angles 
(<Px ,  <p2, y>3)  (Fig.  2).  The  angular  velocity  of  cylinder  is  calculated 
by 

d(px  d(p2  d(p2 

to,  = - ,  to-,  =  — a),  =  — 

dt  ~  dt  dt 


(V) 


and 
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difi 2  dtp2  dfo  d<p2 

W,=tD,,  -  =  -  =  01-,,  -  =F - . 

1  dt  dt  2  dt  dt 


(8) 


hR(fa  +  ^fa,  fa  +  &  fa)  =  MR(fa,  fa) 
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If  (wj  ,oi2,  oi3)  are  given,  time  integration  of  (7)  with  the  time 
interval  At  leads  to 

A  ipl  =  oilAt,  A<p2  =  oiiAt,  A  <p3  =  oi3At.  (9) 

The  increments  (Aif/2,  Ai//2)  are  determined  by  the  relationship  be¬ 
tween  the  two  rotation  matrices  MR(fa  +  Afa,fa  +  Afa)  and 

MR(fa,fa)’ 
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it  is  two-dimensional  with  the  rotation  matrix  given  by 
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2.3  F-Coordinate.  The  F-coordinate  is  represented  by 
FF(X,iF,jF,kF)  with  the  origin  X,  unit  vectors  (iF,jF,kF),  and 
coordinates  (jcF,yF,zF).  Let  Vw  be  the  fluid  velocity.  The  fluid-to- 
cylinder  velocity  is  represented  by  Vr=  Vw-V  that  is  decom¬ 
posed  into  two  parts, 
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Let  the  cylinder  rotate  around  (iF,jF,kF)  with  the  angular  ve¬ 
locity  components  represented  by  ((oF,<nF  ,  roF)  (Fig.  2).  They  are 
connected  to  the  angular  velocity  components  in  the  M-coordinate 
system  by 


(18) 


3  Dynamics 

3.1  Momentum  Balance.  The  translation  velocity  of  the  cyl¬ 
inder  (V)  is  governed  by  the  momentum  equation  in  the 
E-coordinate  system, 
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(19a) 


V,  =  (Vr- iF)ip 

is  the  component  parallel  to  the  cylinder’s  main  axis  (i.e.,  along 
iM),  and 

V2  =  Vr-(V,-iF)iF 

is  the  component  perpendicular  to  the  cylinder’s  main-axial  direc¬ 
tion.  The  unit  vectors  for  the  F-coordinate  are  defined  by  (column 
vectors) 


V2/|V2|,  kF-iFXjF.  (12) 


The  F-coordinate  system  is  solely  determined  by  orientation  of  the 
cylinder’s  main-axis  (iM)  and  the  water-to-cylinder  velocity.  Note 
that  the  M-  and  F-coordinate  systems  have  one  common  unit  vec¬ 
tor  iM  (orientation  of  the  cylinder). 

Let  FR  be  the  rotation  matrix  from  the  F-coordinate  to  the 
E-coordinate, 


4>mf—  OmJf)’  (13) 


(14) 


Here,  i^p  is  the  angle  between  the  two  unit  vectors  (jM  ,jF)-  Let 
the  vector  P  be  represented  by  FP  in  the  F-coordinate.  Transfor¬ 
mation  of  FP  into  FP  contains  rotation  and  translation. 


where  g  is  the  gravitational  acceleration.  II  is  the  cylinder  vol¬ 
ume,  p  is  the  rigid  body  density,  pW=m,  is  the  cylinder  mass,  Fnh 
is  the  nonhydrodynamic  force,  and  F/,  is  the  hydrodynamic  force 
(i.e.,  surface  force  including  drag,  lift,  impact  forces).  Both  F„/, 
and  F/,  are  integrated  for  the  cylinder.  The  drag  and  lift  forces  are 
calculated  using  the  drag  and  lift  laws  with  the  given  water-to- 
cylinder  velocity  (Vr).  In  the  F-coordinate,  Vr  is  decomposed  into 
along-cylinder  (V[)  and  across-cylinder  (V2)  components. 

The  nonhydrodynamic  force  F„/,  is  the  buoyancy  force  (Ffi)  for 
the  air  and  water  phases, 

^nh  =  Fb  =  k(Pang,pwUg), 

where  (pa,Pvl.)  are  the  air  and  water  densities.  F,,;,  is  the  resultant 
of  the  buoyancy  force  (F;,),  pore  water  pressure  force  (F^),  and 
shearing  resistance  force  (Fs)  for  the  sediment  phase  (see  Sec.  6). 

3.2  Moment  of  Momentum  Equation.  It  is  convenient  to 
write  the  moment  of  momentum  equation, 

do.) 

3  —  =  -2J-(DX  m)+Mrf  +  M,„  (20) 

dt 

in  the  M-coordinate  system  with  the  cylinder’s  angular  velocity 
components  (&>j ,  &>2 ,  <u3)  defined  by  (19a)  and  (19Z>).  Here,  the 
first  term  on  the  right-hand  side  is  an  apparent  torque  (similar  to 
the  Corilois  term  in  earth  science)  due  to  the  use  of  the  rotating 
coordinate  system  (i.e.,  the  M-coordinate),  and 

^  =  w23m  +  (21) 

is  the  angular  velocity  of  the  M-coordinate  system.  If  0^=0,  then 
12  =  to,  which  leads  to 


-  2  J  •  (12  X  to)  = 


|  0,  if  toj  =  0  (i.e. ,12  =  to), 

^ —  2,/ot*q +  273to1to2kM,  if  toj  ^  0. 


(15) 


Use  of  the  F-coordinate  system  simplifies  the  calculations  for  the 
lift  and  drag  forces  and  torques  acting  on  the  cylinder.  Since  the 
M-  and  F-coordinates  share  a  common  axis  iM=iF>  the  rotation 
matrix  from  the  F-  to  M-  coordinate  systems  is  given  by 


(22) 

In  this  study,  the  apparent  torque  is  neglected.  The  gravity  force, 
passing  the  COM,  does  not  induce  the  moment.  M„ft  and  MA  are 
the  nonhydrodynamic  and  hydrodynamic  force  torques.  In  the 
M-coordinate  system,  the  moment  of  gyration  tensor  for  the  axi¬ 
ally  symmetric  cylinder  is  a  diagonal  matrix 
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Fig.  3  Three  patterns  of  cylinder  penetration  with  the  cross 
section  being  (a)  a  complete  ellipse,  (b)  a  cutoff  ellipse  with 
one  side  straight  line,  and  (c)  a  cutoff  ellipse  with  two  side 
straight  lines 
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where  Jl,  J2,  and  J3  are  the  moments  of  inertia.  The  buoyancy 
force  induces  the  moment  in  the  jM  direction  if  the  COM  does  not 
coincide  with  the  COV  (i.e.,  xAO), 


M&=|Fi|AfCOS02jM. 


(24) 


Computation  of  nonhydrodynamic  and  hydrodynamic  forces 
(Fnh,Fh)  and  torques  (Mnh,M  h)  is  more  complicated  for  a  cylin¬ 
der  penetrating  through  air-water  and  water-sediment  interfaces 
than  falling  through  a  single  medium  such  as  water.  At  the  in¬ 
stance  when  the  cylinder  penetrates  into  an  interface,  three  situa¬ 
tions  may  exist:  the  cross  section  is  a  complete  ellipse  [Fig.  3(a)], 
a  cutoff  ellipse  with  one  side  straight  line  [Fig.  3(/>)],  or  a  cutoff 
ellipse  with  two  straight  lines  [Fig.  3(c)].  The  interface  separates 
the  cylinder  to  two  parts.  Each  part  contains  a  noncylinder  D  and 
a  subcylinder  C  (Fig.  4).  Let  ( Lc,Ld ),  (Elc,[ld)  and  (nc,II^)  be 
the  lengths,  surfaces,  and  volumes  of  [C,D],  and  (hl,h2)  the 
depths  of  the  two  sides  of  D  (Fig.  5).  The  characteristics  of  the 
geometric  parameters  (Lc,hl,h2)  are  listed  in  Table  1.  The  COV 
for  the  portion  [C,D]  is  called  the  partial  COV  (PCOV). 


4  Equivalent  Cylinder  Method 

4.1  Equivalent  Cylinder.  During  penetration,  the  part  that 
contacts  fluid  (air  or  water)  is  treated  as  a  cylinder  [E]  with  the 
same  mass  and  PCOV  location  and  with  the  assumption  that  the 
buoyancy  and  hydrodynamic  forces  and  torques  for  [C,D]  are  the 
same  for  [E],  The  cylinder  [E],  called  the  equivalent  cylinder,  is 
used  to  represent  the  part  [C,D],  Thus,  the  theoretical  procedure 
developed  for  calculating  external  forcing  (buoyancy  and  hydro- 
dynamic  forces  and  torques)  for  a  cylinder  [8]  can  be  easily  used 
for  [E], 


Fig.  4  Illustration  of  PCOV  (B~),  xu  and  C  for  the  tail  part 
[tfF,#1*]  for  the  case  in  Fig.  3(a) 


Fig.  5  Geometry  of  the  part  D(1) 


4.2  Volume  of  [C,D\.  In  the  M-coordinate  system,  the  area 
of  the  vertical  cross  section  of  D  is  given  by 


h(x) 


j(jc)  =  R~  cos  Ml- - ]-[/?-  h(x)\\lR~  -  [1?  -  h(x)]~. 


(25  a) 


where  h(x)  is  the  depth  of  the  cross  section, 

A/7 

h(x)  =  lil  +  — (x-xY,  A/7  —  h2  —  // 1 , 

Lj 

where  Ld  is  the  length  of  D  (see  Fig.  5).  Integration  of  j(x)  along 
x  axis  gives  the  volume  of  D , 


(25/7) 


If, 


f 

J  x. 


R^L 

s(x)dx  =  =  T rR2ld,  (26) 

A« 


where 


h ,  h2 

k,  =  1-— ,  k2  =  l  —  (27a) 

1  R  2  R 

P(k1,K2)  =  COS^M^l)  -  V  I  -  *1  +  |(1  -  /Cj)3/2  -  K2  COS _1(/f2) 

+  Vl-/c|-|(l  -  a:1)3/2,  {21b) 

and 

RL 

li=^K^-  (27c) 

Here,  ld  is  the  equivalent  length  of  D.  The  volume  of  C  is  calcu¬ 
lated  by 

nc=T tR2Lc.  (28) 

The  total  volume  of  [C,D]  is 

II  =  1 tR2/, 

and 


/  -Lc  +  ld 

is  the  length  of  the  equivalent  cylinder  E. 


Table  1  Geometric  parameters  during  the  cylinder  penetration 


4 

hi 

hi 

Upper  and  lower  parts  of  Fig.  3(a) 

>0 

2R 

0 

Upper  part  of  Fig.  3(b) 

>0 

2R 

0~2R 

Lower  part  of  Fig.  3(b) 

0 

0~2R 

0 

Upper  and  lower  parts  of  Fig.  3(c) 

0 

0~2R 

0~2R 
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4.3  PCOV  of  [C,D\.  Let  (£+,  rf  )  and  (£“,  rf)  be  the  PCOV 
of  the  head  [C,D]  (in  the  direction  of  iM)  or  tail  (in  the  opposite 
direction  of  iM)  [C,D]  (denoted  by  B± ,  positive  sign  for  the  head 
part)  in  the  M-coordinate  system, 

1 


(£*.  A)  =  ^[/JT nc(x,z)dv  +  $  $  hij,x,z)dv\ 


1 


ii,  +  n„ 


x,  ±—  ,0 


nc+  [ 

J  X, 


(x,z)s(x)dx 


t  =  X  1 


RAh 


/3(K„/c2)Lrf(l  +  7rAhLc/TlLdl) 
Ah)  1  -  2\  R 


7f=±  sign(cos  'Ll) 


R 


6/3(Kl  ,k2)(1  +  iTAhL^L-/)^1’^’ 


where 


/JLx(ki,K2)  =  “[(2/Cj  -  l)cOS  1 K2  -  (2k]  -  1)C0S  *K]  +  /C[ V 1  -  K| 

1. 


■  ^2\  1  -  *2]  +  ~ [«^2 V(1  -  kIY  ~  /Cl  V(1  -  /Cf)3] 


-(k2V  1  -  ac2  -  Aq  V 1  -  #cj  +  sin  1k2  -sin  *1^) 

8 

/  - 1  - 1  n  2  2\ 

K\\K2  COS  K2  -  K\  COS  /q  +  V 1  —  /q  -  V 1  —  k2) 


/X-(ki,K2)  =  K,  V(  l  -  K|)3  -  /c2\/(l  -  /r2)3  +  ^(/Cl  Vl  - 


5.1  Buoyancy  Force  and  Torque.  The  buoyancy  force  F6  is 
the  product  of  the  air  (or  water)  density  and  volume, 


Fb  =  p(nc  +  I  l,/)k  =  pirR2(Lc  +  Ij)  k . 


Mfc  -  rPCOV  X  Ffc.  (35) 

Substitution  of  (33)  and  (34)  into  (35)  leads  to 

Mi,  =  -  pttR2{Lc  +  ld)  (f  cos  i//2  +  r)  sin  </f2)jM-  (36) 

5.2  Drag  and  Lift  Forces.  The  drag  and  lift  forces  exerted  on 
the  cylinder  is  represented  by 


F/i  -  (^rfi'F  +  Finif  +  F^kp)  +  F/, 


(37) 


(29) 

where  xx  is  defined  as  the  location  of  interface  between  C  and  D. 
Substitution  of  (25 a),  (25 b),  (26),  and  (28)  into  (29)  leads  to 


(30) 


where  {Fdl,Fd2,F A  are  the  components  of  drag  force  along  iF 
(along-cylinder),  jF  (across-cylinder)  and  kF  directions.  F;  repre¬ 
sented  the  lift  force.  Linearization  of  drag  and  lift  laws  is  used  in 
the  computation. 

Let  (Cdl ,  Cd2)  be  the  drag  coefficients  in  the  along-  and  across- 
cylinder  directions  (Reynolds  number  dependent).  The  drag  force 
coefficients  are  calculated  on  the  base  of  steady  flow;  they  is 
different  from  the  fluid  around  an  accelerated  solid  body.  The 
added  mass  correction  is  represented  by  the  ratios  (fl  ,f2,f2)  in  the 
three  directions  of  the  F-coordinate  system. 

The  drag  force  along  iF  is  calculated  by 


(31) 


Fd  1  -  Ctd  1  (?)  V 1 , 


cldM  —  cd 


ttR1 


2  (l+/i) 


|Ft(0l- 


(38) 

(39) 


Cd  1  is  almost  independent  on  the  axial  Reynolds  number  (Re) 
when  Re>  104,  but  dependent  on  the  cylinder’s  aspect  ratio  [9], 


1.0,  if  S>  8, 

0.75  +  <5/32.1934  +  0.09612/ <52,  if  8  &  <5>  0.5, 
1.15,  if  <5=s  0.5. 


Substitution  of  (11)  and  (12)  into  (38)  leads  to 


-  k2 yl  -  k2  +  sin  '/<j-sin  1k2).  (32) 

Let  (f ,  97)  be  (i;+,y+)  for  the  head  part  and  for  the  tail 

part.  The  position  vector  of  PCOV  in  the  M-coordinate  system  is 
represented  by 

rPCOV  =  +  V^-m-  (33) 


5  Forces  and  Torques  in  Air  and  Water 

Calculation  of  the  buoyancy  force  and  torque  is  straightfor¬ 
ward.  Calculation  of  the  surface  force  and  torque  is  not  simple. 
Assume  that  the  surface  force  and  torque  on  the  equivalent  cylin¬ 
der  E  are  the  same  on  the  \C,D],  If  [C,D]  moves  in  fluid  (air  or 
water),  the  recursive  model  recently  developed  [8]  can  be  used  to 
calculate  for  equivalent  cylinder  E.  Thus,  the  water  column  is 
taken  as  the  example  to  illustrate  the  calculation  of  the  hydrody¬ 
namic  force  and  torque.  Computation  of  the  surface  force  and 
torque  due  to  sediment  is  described  in  Sec.  6. 


Fd dF  -  -  Ctdi(f)In  ■ 


1 

u 

\ 

V 

- 

vw 

\ 

w 

/ 

In  -  Mf’ 


(40) 


(41) 


where  the  superscript  “T”  denotes  the  transpose. 
The  drag  force  along  jF  is  calculated  by 


4 


L/2-x 

U2-X 


CAVtfj^-^dV=  Ctd2(t)V2  +frd2(t). 


(42) 


where 


V2(v)  =  V2-w2rj 

is  the  water-to-cylinder  velocity  at  the  surface  in  the  jF  direction 
and 


Cld2(t)  —  2  Cd2LR 


(1  +m  2 


V- 


2  F 


(43a) 


The  torque  due  to  the  buoyancy  force  for  the  upper  or  lower  part 
is  given  by 
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(34)  frAt)  -  2 C^LR^jzjy-A  +  ^2)(4)2-  (43*) 


An  empirical  formula  is  used  for  calculating  Cd2  [10] 
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c,t 


1.9276  +  8/Re, 
1.261  +  16/Re, 
0.855  +  89/Re, 

0.84  +  0.000  03  Re, 
1.2  -  4/(5, 

0.835  -  0. 35/5, 

0.7  -  0.08/ <5, 
1.875-0.000  004  5  Re, 
1/(64 1, 550/Re  +  1.5), 


if  Re  12, 
if  12  <  Re  ^  180, 
if  180  <  Re  =£  2000, 
if  2000  <  Re  «  12,000, 
if  12,000  <  Re  «  150,000,  <53=  10, 

if  12,000  <  Re  =S  150,000,  2  *£  <5<  10, 

if  12,000  <  Re  «  150,000,  <5<  2, 

if  150,000  <  Re  «  350,000, 
if  Re  >  350,000. 


(44) 


Substitution  of  (11)  and  (12)  into  (42)  leads  to 


Fcni f  -  -  Ctd2(t)I22  ' 


/ 

u 

\ 

V 

- 

Vw 

\ 

w 

Ww 

/ 

+  /<-rf2(0jF i  122-jpjF- 


velocity  in  the  kF  direction, 

¥  /  F 

^3  =  W2  V- 

The  drag  force  along  kF  is  calculated  by 

-L 

°2 1  ^2  * 


3  ■ 


where 


s~>  r>  Pw  F\  F\ 
Cd2^~ - —0)i\(O' 


(i  +/2) 

frdM^F, 


j2~*r?dV-j^  rf  dr\ 


frd 30  -  -  ^Crf2 77^77 x(3i2  +  4^2)|(uf|u)f 

6  (l+/2) 


is  the  rotational  drag  force  in  the  kF  direction. 

The  water-to-cylinder  velocity  determines  the  lift  force  [11] 


C„(t) 


rL/2-x 

I  1 

J  -L/2~x 


V'2(y)dy 


kF,  CM)  —  CiLR—^—\V2\. 

U  +J2) 


■CM)  i 


32  ' 


/ 

u 

uw 

\ 

V 

- 

Vw 

\ 

w 

Ww 

/ 

+/w(r)kF,  I32  —  kFjF, 


where 


fM)  —  CM)X“  3 


F;,  -  -  [Ct,i(/)In  +  C,d2{t) l22  +  C,,(f)I32] 


1 

u 

uw 

\ 

V 

- 

Vw 

\ 

w 

Ww_ 

/ 

(45) 

The  angular  velocity  (<u2)  causes  nonuniform  water-to-cylinder 


r12 

/ 

r13 

X 

r22 

+  LfrdM)  +frl(t)l 

r23 

_r32_ 

r 

/a  _ 

+/ne(f) 


(52) 


Substitution  of  (52)  into  (19a)  leads  to  the  cylinder’s  momentum 
equation, 


(46) 


D 


+  «i> 


(19*) 


where 


(47) 


(48) 


0 

r 

r\2 

r'n 

uw 

- 

0 

+  Z?1 

r22 

+  b 9 

r23 

Ww_ 

(1  -pjp)g 

_r32_ 

_r33  _ 

=  D 


r,  Ctdl(t)ln  +  CldM)I22  +  CM)h3  ,  frdl^f) 

D  = - — - ,  *[  = 


PU 


pU 


bn  — 


frilM)  +fM) 

Pn 


(49) 


where  C;  is  the  lift  coefficient.  An  empirical  formula  is  used  for 
calculating  C;  [12], 

[  2co,R/V2,  if(u,R/R2«4, 

C  =  1  1  2  (50) 

1  [s  +  O^f&qR/Ro-T),  if  co1R/V2  >  4. 

Substitution  of  (11)  and  (12)  into  (49)  leads  to 


5.3  Drag  and  Lift  Torques.  For  an  axially  symmetric  cylin¬ 
der,  the  moment  of  the  hydrodynamic  force  in  the  iF  direction  is 
not  caused  by  the  drag  and  lift  forces,  but  by  the  viscous  fluid.  The 
moment  of  the  viscous  force  of  steady  flow  between  two  rotating 
cylinders  with  the  common  axis  is  calculated  by  [1] 


M 


rx 

-  4 17 fi  — 


l("i  ~"o), 


(51) 


where  (rx,r0)  and  (&>j ,  co0)  are  the  radii  and  angle  velocities  of  the 
inner  and  outer  cylinders  and  p  is  the  viscosity.  The  moment  of 
the  viscous  force  on  one  rotating  cylinder  is  th  limit  case  of  the 
two  rotating  cylinders  as  r0—>°o  and  co0=0.  The  moment  of  the 
viscous  force  around  iF  is  calculated  by 

M«i  =  -  C^itoiip,  Cm  1  =  i Tfj,Ld2-  (53) 

Same  as  the  hydrodynamic  forces,  the  torques  along  the  jF  and 
kF  axes,  (Mdl  ,Md2,Mf),  are  calculated.  When  the  cylinder  rotates 
around  jF  with  the  angular  velocity  co^,  the  drag  force  causes  a 
torque  on  the  cylinder  in  the  jF  direction, 

rLI2-x 


Mj 


is  the  rotational  lift  force.  Substitution  of  (41),  (45),  (47),  and  (51) 
into  (37)  and  use  of  (14)  lead  to 


Fl  F| 
■  a)-,\(Oi\ 


[ 


C,nR 


L/2-X 


(1  +/,) 


y2\y\dy 


Jf 


:  _  [C„,2(t)&>2]jF, 


(54) 
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Pw 


Cm2{t)  —  9  CdlR 


i  — L4  +  —L2^  +  X* 

(1  +/r)V  16  2 


141 


where  fr  is  the  added  mass  factor  for  the  moment  of  drag  and  lift 
forces.  If  the  water-to-cylinder  velocity  or  the  cylinder  mass  dis¬ 
tribution  is  nonuniform  (^AO),  the  drag  force  causes  a  torque  on 
the  cylinder  in  the  kF  direction, 
rLI2-x 


Mj 


l 


C,nR 


-U2-X 


(1  +fr) 


(V2-  io3y)-ydT] 


'■  ~  [Cm3(^)4  +  ^3(f)]kF> 


(55) 


Cm3(t)  -  CcP.Rj~~rjy  \  V2L3  +  V2L)c  +  l-L3<oF3x  +  L^J , 


(56  a) 


M3(t)  =  Cd2R7^~V22LX. 


(56  b) 


(1  +/,) 

The  lift  force  exerts  a  torque  on  the  cylinder  in  the  jF  direction, 

I'M-* 

C/R~(V2  -  W37])T}dT) 

J kr 


Mr 


r 

j  -i 


L/2-x 

\Crn[(t)(i>3  +  Mi(t)\ jF, 

1 


Jf 


(57) 


(l+/r)  V  12 


Cjt)  -  ClV2R7f^-L\  i-L2  +  r),  M,(t)  =  R^LVlx- 


f kr 


(58) 

After  the  angular  velocity  components  (coF,  (o3)  are  transformed 
into  (i o2,co3 )  (from  the  F-coordinate  to  the  M-coordinate)  using 
(18),  and  the  unit  vectors  (jF,kF)  are  transformed  into  (jM,kM) 
using  the  rotation  matrix  (17),  the  drag  force  torques  in  the  jF 
direction  (54)  and  in  the  kF  direction  (55)  are  represented  by 


Md2=-cm2m. 


22  ' 


(o2 

L"3  J 


,  H22  =  e2e2, 


(59) 


-  -  Cm3(/)H33  • 


(02 

L^3  J 


-  M3(t)e3,  H33  =  e3e3 ,  (60) 


and  the  lift  torque  in  the  jF  (57)  is  represented  by 


M/2  -  C„,/(r)H23  ■ 


+  M/(r)e2,  H23-e2e3.  (61) 


co2 
0)3 

Summation  of  (53)  and  (59) — (6 1 )  leads  to 
M;,  =  M„  +  +  Mf/3  +  M/2 

=  _  Cm|Ui|iF -  [C,„2(/)H22  +  Cm3(r)H33  -  Cm;(r)H23]  ■ 

+  M,(f)e2-M3(r)e3. 


(o2 

L"l  J 
(62) 


6  Resistant  Forces  in  Sediment 

6.1  Water  Cavity.  As  the  cylinder  impacts  and  penetrates 
into  the  sediment,  it  pushes  the  sediment  and  leaves  space  in  the 
wake.  This  space  is  refilled  by  water  right  away  and  a  water  cavity 
is  produced  (Fig.  6).  At  the  instant  of  the  penetration,  the  total 
resistant  force  on  the  cylinder  is  represented  by 

F  =  J  VS{f°b  +  iJ  +  rb+f£\dar+¥pw,  (63) 

where  (f^.f,;,)  and  (fJJ  ,fj,1)  are  the  sediment  buoyancy  and  shear 
resistance  forces  and  water  buoyancy  and  hydrodynamic  forces 


Fig.  6  The  impact  (resistant)  force  exerted  on  the  part  of  the 
object’s  surface  moving  towards  the  sediment 


(per  unit  area)  at  the  point  r  over  the  cylinder’s  surface;  ased  is  the 
area  of  the  cylinder’s  surface  below  the  water-sediment  interface; 
and  Fpvi,  is  the  pore  water  pressure  force  on  the  whole  cylinder.  In 
the  sediment,  the  magnitude  of  the  sediment  nonhydrostatic  force 
is  much  larger  than  the  magnitude  of  the  water  hydrodynamic 
force, 

|f'1  >  |f£|, 

which  means  that  f)|  in  (63)  can  be  neglected.  The  water  buoyancy 
force  per  unit  area  over  the  cylinder’s  surface  is  defined  by 

fr  =  -p»f?4s-z)n,  (64) 

where  zws  is  the  depth  of  the  water-sediment  interface  and  n  is  the 
unit  vector  normal  to  the  cylinder  surface  (outward  positive). 

Let  v  be  the  velocity  at  point  r  (represented  in  the 
M-coordinate)  on  the  cylinder  surface, 

v  =  V  +  a)  X  r. 


The  step  function  S  is  defined  by 


f  1,  v  •  n  >  0, 
[o,  v  •  n  =S  0, 


(65) 


which  shows  that  the  sediment  buoyancy  and  shear  resistance 
forces  act  when  the  cylinder  moves  towards  it.  Let  v„  be  the 
normal  velocity.  The  tangential  velocity  is  represented  by 


VT=V-V„.  (66) 

The  tangential  unit  vector  (r)  is  defined  by 


vT 


which  is  opposite  to  vr  (Fig.  7). 


(67) 


6.2  Sediment  Resistant  Forces.  When  the  cylinder  impacts 
and  penetrates  into  the  sediment,  it  will  create  a  large  transient 
pore  pressure  in  the  sediment  that  causes  ruptures  in  the  sediment 
which  influences  the  lifting  forces  on  the  cylinder  [13,14], 

The  sediment  buoyancy  force  per  unit  area  is  defined  by 


pM)gdz' , 


(68) 


where  ps{z )  is  the  sediment  density. 

The  shear  resistant  force  per  unit  area  fs;,  depends  on  the  cyl¬ 
inder’s  penetration  speed  ( V )  and  the  sediment  strength.  Let  S(z ) 
be  the  sediment  shear  strength.  The  shear  strength  is  defined  as  the 
maximum  stress  that  a  material  can  withstand  before  failure  in 
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♦  V 


Fig.  7  Momentum  and  angular  momentum  balance  for  the  cyl¬ 
inder’s  penetration  through  the  water-sediment  interface 


shear.  Calculation  of  shear  strength  depends  upon  the  test  method. 

After  entering  the  water-sediment  interface,  the  cylinder  re¬ 
duces  its  speed  (V),  and  the  sediment  shearing  resistant  force  also 
decreases.  When  the  cylinder  ceases,  the  shearing  resistant  force 
should  be  the  same  as  the  sediment  shear  strength  S(z ).  Thus,  the 
shearing  resistant  force  is  represented  by 


tsh  =  S(z)G(V)T,  G(0)=1, 
where  G(V)  is  the  impact  function  defined  by 


G(V)=A 


1  -(1  -A_1)expl 


(69) 

(70) 


Here,  Vrest  is  an  infinitesimally  small  value  for  V  representing  the 
cease  of  the  cylinder  in  the  sediment.  The  impact  function  has  the 
following  feature, 


Lim  G(V)  =  A,  (71) 

v— >°° 


which  shows  that  when  the  cylinder  impacts  on  the  sediment  (usu¬ 
ally  with  large  penetration  speed),  the  impact  function  takes  the 
value  of  A.  Thus,  we  may  call  A  the  impact  factor.  Note  that  A  and 
Vfest  are  the  two  tuning  parameters  of  the  numerical  model.  In  this 
study  we  use 

A  =  10,  Vrest  =  0.04  m  s-1.  (72) 


The  shear  strength  of  the  sediment  is  directly  measured  from  the 
gravity  cores  using  the  fall  cone  apparatus  (model  G-200)  (see 
Sec.  8.2). 

The  total  force  due  to  the  pore  water  pressure  on  the  cylinder  is 
computed  by  [15] 


F„ 


(T^  + 


1  +  e„  dw  \ 


- 


dt  / 


(73) 


kp  is  the  permeability  coefficient  (10' 


'  4  m  s  1  [15]), 


where 

e„(~0.50)  is  the  void  ratio,  and  B  is  the  length  of  the  rupture  line. 
Substitution  of  (64),  (68),  (69),  and  (73)  into  (63)  leads  to 


1 


t {SfjiG(V)S(z)]d, 


r-J  »[(*[ 


ps{z')gdz‘ 


+  Pwg(zws  ~  Z) 


77 

dcr+k-ps(z) 

O 


1  +  e„  dw 
c v  dt 


B3, 


(74) 


which  is  the  external  force  acted  on  the  cylinder  in  the  sediment 
phase.  The  torque  due  to  the  sediment  (Ms)  is  calculated  by 


M 


■L 


(r  X  T)[SfjLG{V)S{zY\dcr+  (r  X  n) 


1 


8  Ps(z')gdz'  +  pwg(zm  -  z) 


da 


77  /  gW 

+  (rpw  X  k)-ps(z)(  —  + 


1  +  e„  dw 
ev  dt 


B3. 


(75) 


where  rpw  is  the  position  vector  (in  the  M-coordinate)  indicating 
the  location  of  the  cylinder’s  rupture  line. 


7  Model  Integration 

The  momentum  equation  (19a)  and  (197;)  and  moment  of  mo¬ 
mentum  equation  (20)  are  integrated  numerically  using  the  triple 
coordinate  transformation.  The  momentum  equation  is  integrated 
in  the  E-coordinate  system.  The  hydrodynamic  (drag  and  lift) 
force  is  transformed  from  the  F-coordinate  to  the  E-coordinate. 
The  moment  of  momentum  equation  is  integrated  in  the 
M-coordinate  system.  The  hydrodynamic  torque  is  transformed 
from  the  F-coordinate  to  the  M-coordinate.  After  the  cylinder  pen¬ 
etrates  into  the  sediment,  the  resistant  force  due  to  sediment  f1 
reduces  the  cylinder’s  speed  and  changes  the  turning  angle. 


7.1  Cylinder’s  Angular  Velocity.  Substitution  of  (24)  and 
(62)  into  (20)  leads  to  the  equations  for  (wl,a)2,(o 3). 


where 


B  = 


d(o{ 

dt 


(Oi 


dt  _ 


■«t"t, 


(o2 

aij 


+  a2. 


(76) 

(77) 


Cm  1  81 rfj.L 


J , 


Pn 


-  0 


0  - 


'  LCm2(t)H22  +  Cm3(?)H33  -  Cm/(i)H23], 


-  0 


0  - 


(M,e2  -  M3e3)  +  ^y^cos  i/i2 


■  (78) 


Equation  (76)  has  the  analytical  solution 

&)I(f)  =  &>i(i0)exp[-a1(f®70)],  (79) 

which  represents  damping  rotation  of  the  cylinder  around  the 
main  axis  (iM).  Substitution  of  (79)  into  (8)  leads  to 


dtp  1 

—  =  oq^expt-  a.\(t—  f0)], 
dt 


and  its  integration  leads  to 
co .  (7n) 

<Pi(f)  =  - - -exp[-  ai(f-  f0)l  +  <Pt(fo)-  (80) 

Cl  J 

Equations  (79)  and  (80)  are  the  analytic  formulas  for  predicting 
the  angle  and  angular  velocity  around  the  cylinder’s  main  axis 
(“i.¥>t)- 


7.2  Recursive  Procedure.  The  basic  equations  (19a),  (19 b), 
(77),  (79),  and  (80)  describe  the  dynamics  of  the  falling  cylinder. 
It  is  noted  that  the  coefficient  matrices  B,  D  and  the  vectors  ax. 
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a2  depend  on  drag/lift  coefficients.  Besides,  B,  D,  or,,  a2  depend 
on  the  fluid-to-cylinder  velocity  and  cylinder’s  angular  velocity. 
Equations  (19a),  (19 b),  and  (79)  are  nonlinear  equations. 

Let  matrices  B  and  D  be  separated  into  diagonal  and  nondiago¬ 
nal  parts, 


1 

0 

O 

o 

d\2 

d» 

0 

d2 

0 

III 

Q 

d2l 

0 

dn 

1 

O 

0 

m 

_d2i 

dyi 

1 

o 

(81) 


bx 

0  " 

0  bn 

Bi  = 

,  b2  = 

.0 

b2. 

b2\  0  _ 

Substitution  of  (81)  into  (19a)  and  (19 b)  leads  to 


(82) 


—  =  -  V  +  /?,  p  =  ai  -  D2  •  V, 


U 

Vi 

V  = 

V 

= 

v2 

w 

Vl- 

(83a) 


and  substitution  of  (82)  into  (77)  leads  to 


fc>2 

=  -B, 

<x>2 

+  % 

1 

<§’ 

III 

(L>2 

."3. 

6l>3 

co3 

(83b) 


IfB,,  Db  p,  y  are  taken  the  values  given  at  the  present  time  step 
f,„  (83a)  and  (83Z>)  can  be  treated  as  “linear”  equations  (local 
linearization)  on  [f„,f„  +  Af]  and  integrated  analytically, 


and 


V&n+ 1)  ' 


Vi(tn)  +  PAQAt,  if  dj(tn)  =  0, 

if  di(tn)  #  0,  i=  1,2,3, 

",(0  +  Ji(tn)At,  if  dj(t„)  =  0 , 

",(0  -  rzr) exp[-  HOAt]  +  7777 ,  if  <7,(0  ^  0,  1  =  2,3. 


Integration  of  (83a)  twice  from  tn  leads  to  the  translation  position  of  the  cylinder  at  fn+1, 

Xi(tn)  +  Vi(t„)At  +  \_Pi(t„)Ar,  if  di(tn)  =  0, 


Xi(tn+ 1) : 


,  x  A (O.  1 

Xi(t„)  +  —At  -  ■ 


Vi(t„ ) 


AW 


l{exp[-  di(tn)At]  -  1},  if  di(t„)  +  0, 


di(tn)  dj(t„)  \  '  "  di(tn), 
where  x2=x,  x2  =  y,  and  .r3  =  z.  Integration  of  (83 b)  twice  from  t„  leads  to  the  change  of  rotation  angles  (<p2,  <p3)  at  tn+l, 

Mi(tn)At+  \yi(tn)At2,  if  bi(tn)  =  0, 


A<Pi(tn+ 1)  = 


“i(0 


■Ar 


.  bi(tn)  bi(tn) 


Viitn)  ~ 


bi(t„ ) 


l{exp[-  b;(tn)At]  -  1},  if  b,(t„)  ^0,  i  =  2,3 . 


(84) 


(85) 


(86) 


(87) 


Equations  (84)  and  (85)  are  the  recursive  formulas  for  predict¬ 
ing  the  cylinder’s  translation  velocity  ( u,v ,w )  and  angular  veloc¬ 
ity  (iol,w2,io2),  and  Eqs.  (86)  and  (87)  are  the  recursive  formulas 
for  predicting  the  cylinder’s  translation  position  (x,y,z)  and  rota¬ 
tion  angle  increments  ( A<p2,Aip3 )  in  the  M-coordinate  system.  The 
cylinder’s  orientation  is  represented  by  angles  (<//| ,  ifi2,  ip2)  with 
i/f|  =  (P[,  and  a  relationship  between  (Ai//2,Ai//3)  and  (A<p2,A<p3) 
given  by  (10). 

Let  [x(t0)  ,y(t0)  ,z(to) ,u(t0) ,u(to)  ,w(t0)]  be  the  cylinder’s  initial 
translation  and  velocity  and  ,i[/2(t0)  ,i//2(t0)  ,a>i(t0) , 

a)2(f0),a)3(r0)]  be  the  cylinder’s  initial  orientation  and  angular  ve¬ 
locity.  Following  the  procedures  listed  in  Fig.  8,  the  values  of 
these  variables  for  next  time  step  (f=/q)  are  calculated.  Repeat  of 
the  procedures  leads  to  predicting  the  cylinder’s  position  and  ori¬ 
entation  as  falling  through  the  water  column. 

Theoretically,  the  model  integration  stops  when  the  vertical  co¬ 
ordinate  of  COM  [i.e.,  z(f)]  in  the  E-coordinate  and  the  elevation 
angle  if/2(t)  in  the  M-coordinate  do  not  change  with  time  (in  the 
sediment  column): 


(88) 


Practically,  the  following  criteria  are  used  to  stop  the  integration, 


dz 

<#2 

dt 

55  £l, 

dt 

where  (e!,e2)  user-defined  small  values.  In  this  study,  we  use 
£[  =  l(L6m,  e2  =  10-4. 

8  Cylinder  Drop  Experiments 

Two  cylinder  drop  experiments  were  conducted  to  collect  data 
for  the  model  evaluation.  Exp-1  was  designed  to  collect  data  on  a 
cylinder’s  motion  in  the  water  column  for  various  combinations  of 
the  cylinder’s  parameters.  Exp-2  was  designed  to  collect  synchro¬ 
nized  data  on  sediment  parameters  (shear  strength  and  density) 
and  the  cylinder’s  burial  depth  and  orientation. 

8.1  Exp-1.  Exp-1  was  conducted  at  the  Naval  Postgraduate 
School  swimming  pool  in  June  2001  [16].  It  consisted  of  dropping 
each  of  three  model  cylinders  (Fig.  9)  into  the  water  where  each 
drop  was  recorded  underwater  from  two  viewpoints.  The  physical 
parameters  of  the  model  cylinders  are  listed  in  Table  2.  Figure  10 
depicts  the  overall  setup.  The  controlled  parameters  for  each  drop 
were  LI R  ratio,  ^--value,  initial  velocity  (Vjn),  and  drop  angle.  The 
E-coordinate  system  is  chosen  with  the  origin  at  the  corner  of  the 
swimming  pool  with  the  two  sides  as  x  and  y  axes  and  the  vertical 
z  axis.  The  initial  injection  of  cylinders  was  in  the  (y,z)  plane 
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Fig.  8  Procedure  of  the  recursive  model 


(Fig.  11). 

Initial  velocity  (Vjn)  was  calculated  by  using  the  voltage  return 
of  an  infrared  photo  detector  located  at  the  base  of  the  cylinder 
injector.  The  infrared  sensor  produced  a  square  wave  pulse  when 
no  light  was  detected  due  to  blockage  caused  by  the  cylinder’s 
passage.  The  length  of  the  square  wave  pulse  was  converted  into 
time  by  using  a  universal  counter.  Dividing  the  cylinder’s  length 
by  the  universal  counter’s  time  yielded  Vin.  The  cylinders  were 
dropped  from  several  positions  within  the  injector  mechanism  in 
order  to  produce  a  range  of  Vin.  The  method  used  to  determine  Vin 
required  that  the  infrared  light  sensor  be  located  above  the  water’s 
surface.  This  distance  was  held  fixed  throughout  the  experiment  at 
10  cm. 

The  drop  angle  (initial  value  of  iji"11)  was  controlled  using  the 


Fig.  9  Internal  components  of  the  model  cylinder 
Journal  of  Applied  Mechanics 


Table  2  Physical  parameters  of  the  model  cylinders 


Cylinder 

Mass 

(g) 

L 

(cm) 

Volume 

(cm3) 

Pm 

(g  m-3) 

J 1 

(g  m2) 

X 

(cm) 

^2(^3) 

x(g  m2) 

322.5 

15.20 

191.01 

1.69 

330.5 

0.00 

6087.9 

i 

0.74 

5783.0 

1.48 

6233.8 

2 

254.2 

12.10 

152.05 

1.67 

271.3 

0.06 

3424.6 

0.53 

3206.5 

1.00 

3312.6 

3 

215.3 

9.12 

114.61 

1.88 

235.0 

0.00 

1695.2 

0.29 

1577.5 

0.58 

1556.8 

drop  angle  device.  Five  screw  positions  marked  the  15  deg,  30 
deg,  45  deg,  60  deg,  and  75  deg.  The  drop  angles  were  determined 
from  the  lay  of  the  pool  walkway,  which  was  assumed  to  be  par¬ 
allel  to  the  water’s  surface.  A  range  of  drop  angles  was  chosen  to 
represent  the  various  entry  angles  that  air-  and  surface-laid  mines 
exhibit  in  naval  operation.  This  range  produced  velocities  whose 
horizontal  and  vertical  components  varied  in  magnitude.  This  al¬ 
lowed  for  comparison  of  cylinder  trajectory  sensitivity  with  the 
varying  velocity  components. 

For  each  drop  the  cylinder  was  set  to  a  x  value.  For  positive  \ 
value,  the  cylinders  were  placed  into  the  injector  so  that  the  COM 
was  located  below  the  geometric  center.  For  negative  \  value,  the 


Fig.  10  Exp-1  equipments 
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Table  3  Number  of  drops  conducted  for  different  drop  angles 
and  x  values  for  LI Ft=  15/2 


4i") 

15  deg 

30  deg 

45  deg 

60  deg 

75  deg 

Xi 

13 

15 

15 

15 

12 

Xi 

9 

15 

15 

15 

9 

Xo 

12 

14 

15 

18 

6 

X-i 

0 

6 

6 

6 

0 

X-2 

2 

6 

6 

0 

0 

Table  4  Measuring  ranges  of  the  gravity  cores 


Weight 

(g) 

Apex- 

angle 

Penetration 

(mm) 

Undrained 
shear  strength 
(kPa) 

400 

30  deg 

4.0-15.0 

25-1.8 

100 

30  deg 

5.0-15.0 

4-0.45 

60 

60  deg 

5.0-15.0 

0.6-0.067 

10 

60  deg 

5.0-20.0 

0.10-0.0063 

COM  was  located  above  the  geometric  center  to  release.  A  series 
of  drops  was  then  conducted  in  order  of  decreasing  mine  length 
for  each  angle.  Table  3  indicates  number  of  drops  conducted  for 
different  drop  angles  and  \  value  for  L/R=  15/2.  The  number  of 
drops  for  other  L/R  ratios  (12/2,9/2)  is  comparable  to  that  for 
the  L/R  ratio  of  15/2.  All  together  there  were  712  drops.  Each 
video  camera  had  a  film  time  of  approximately  1  h.  At  the  end  of 
the  day,  the  tapes  were  replayed  in  order  to  determine  clarity  and 
optimum  camera  position. 

Upon  completion  of  the  drop  phase,  the  video  from  each  cam¬ 
era  was  converted  to  digital  format.  The  digital  video  for  each 
view  was  then  analyzed  frame  by  frame  (30  Hz)  in  order  to  de¬ 
termine  the  mine’s  position  in  x-z  and  y-z  planes.  The  cylinder’s 
top  and  bottom  positions  were  input  into  a  MATLAB  generated 
grid,  similar  to  the  ones  within  the  pool.  The  first  point  to  impact 
the  water  was  always  plotted  first.  This  facilitated  tracking  of  the 
initial  entry  point  throughout  the  water  column.  The  cameras  were 
not  time  synchronized;  thus,  the  first  recorded  position  corre¬ 
sponded  to  when  the  full  length  of  the  mine  was  in  view. 

8.2  Exp-2.  Exp-2  was  conducted  on  the  R/V  John  Martin  on 
May  23,  2000  [17].  The  barrel  with  density  ratio  of  1.8  was  re¬ 
leased  horizontally  while  touching  the  surface.  The  initial  condi¬ 
tions  are 

Vin  =  0,  4m)  =  90  0  ,  (90) 

This  would  be  to  eliminate  any  chance  of  inertial  effects  caused 
by  uneven  introduction  into  the  air-sea  interface.  This  also  set  the 
initial  velocity  parameter  in  the  code  to  zero.  The  barrel  was  to  be 
released  17  times.  The  diver  would  snap  the  quick-release  shackle 
on  the  barrel  and  then  dive  down  to  conduct  measurements.  The 
average  depth  of  the  water  was  13  m.  Since  it  was  uncertain  the 
path  the  barrel  would  follow,  both  the  releasing  diver  and  a  second 
safety  diver  would  stay  on  the  surface  until  after  the  barrel  had 
dropped.  Once  reaching  the  bottom,  one  diver  would  take  penetra¬ 
tion  measurements  using  a  meter  stick  marked  at  millimeter  incre¬ 
ments  while  the  other  would  take  a  gravity  core.  After  17  drops, 
the  divers  began  to  run  out  of  air  and  results  were  not  varying 
greatly  so  the  decision  was  made  to  end  the  experiment.  Upon 
return  to  the  Monterey  Bay  Aquarium  Research  Institute,  the  grav¬ 
ity  cores  were  taken  immediately  to  the  USGS  Laboratories  in 
Menlo  Park,  CA  where  they  were  refrigerated  until  the  analysis 
could  be  performed  on  May  31-June  1,  2000. 

Analysis  of  the  gravity  cores  was  begun  on  May  31,  2000  at  the 
USGS  Laboratories  in  Menlo  Park,  CA.  The  gravity  cores  were 
sliced  into  2  cm  segments  to  a  depth  of  10  cms,  and  then  sliced 
into  4  cm  segments.  A  fall  cone  apparatus  (Model  G-200)  was 
used  to  determine  sediment  density  ps(z)  and  shear  strength.  In  the 
test,  it  is  assumed  that  the  shear  strength  of  sediment  at  constant 
penetration  of  a  cone  is  directly  proportional  to  the  weight  of  the 
cone,  and  the  relation  between  undrained  shear  strength  s  and  the 
penetration  h  of  a  cone  of  weight  Q  is  given  by 

S(z)  =  KQ/h2,  (91) 


ing  the  measuring  range  listed  in  Table  4.  The  cones  are  sus¬ 
pended  from  a  permanent  magnet.  By  pressing  a  knob,  the  magnet 
is  moved  so  that  the  magnetic  field  is  broken  momentarily  and  the 
cone  is  released.  Measurements  are  taken  of  penetration  depth  and 
the  evolution  is  repeated  five  times  per  sediment  slice.  These  val¬ 
ues  are  then  averaged  and  correlated  with  a  table  which  gives 
shear  strength.  Previous  studies  [18]  showed  that  the  sediment 
parameters  are  the  most  critical  element  in  determining  how  deep 
an  object  was  buried  when  it  came  to  rest.  During  the  experiment 
at  the  Monterey  Bay,  we  obtained  17  gravity  cores.  Sediment  bulk 
density  and  shear  strength  profiles  (Fig.  12)  generally  show  in¬ 
crease  with  depth  until  approximately  6-9  cm  below  the  water- 
sediment  interface. 

9  Model-Data  Comparison 

The  U.S.  Navy  has  a  2D  model  in  the  x-z  plane  (IMPACT28)  to 
predict  the  cylinder’s  trajectory  and  impact  burial.  Since  the  mo¬ 
tion  of  the  cylinder  is  3D,  the  impact  burial  prediction  using  the 
2D  model  has  large  errors  [19-21],  In  this  study,  a  new  3D  model 
(called  IMPACT35)  is  developed  on  the  base  of  momentum  bal¬ 
ance  (19a)  and  ( 19£>)  and  moment  of  momentum  balance  (20) 
using  a  triple  coordinate  transform  [8]  and  cylinder  decomposi¬ 
tion.  To  evaluate  the  value  added  to  the  3D  model,  comparison 
among  the  observed  data  (from  Exp-1  and  Exp-2)  and  predicted 
data  using  2D  (IMPACT28)  and  3D  (IMPACT35)  models  is  con¬ 
ducted.  Since  position  and  orientation  of  the  cylinder  were  tape 
recorded  after  it  is  submerged  into  the  water,  the  free  water  sur¬ 
face  effect  was  not  detected  from  Exp-1  and  Exp-2. 

9.1  Comparison  Using  Exp-1  Data.  Improvement  from  IM- 
PACT28  to  IMPACT35  in  predicting  cylinders’  trajectory  and  ori¬ 
entation  in  the  water  column  is  verified  using  the  Exp-1  data. 
Here,  we  list  two  examples. 

Positive  x  (Nose-Down):  Cylinder  1  (L=  15.20  cm,  p 

=  1.69  g  cm-3)  with  ^=0.74  cm  is  injected  to  the  water  with  the 


where  K  is  a  constant  which  depends  mainly  on  the  angle  of  the  Fig.  12  Mean  sediment  density  Ps{z)  and  shear  strength  S(z) 
cone,  but  is  also  influenced  by  the  sensitivity  of  the  clay/sediment,  profiles  in  the  Monterey  Bay  collected  during  the  cylinder  drop 
Four  different  cones  are  used  with  this  instrument,  each  one  hav-  experiment  on  May  31,  2000 
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Experiment 


3-D  Model 


2-0  Model 


Fig.  13  Movement  of  cylinder  1  (i=15.20  cm,p=1.69  g  cm-3)  with  *=0.74  m  and  drop  angle  45 
deg  obtained  from  (a)  experiment,  (b)  3D  IMPACT35  model,  and  (c)  2D  Impact28  model 


drop  angle  45  deg.  The  physical  parameters  of  this  cylinder  are 
given  by 

m  =  322.5  g,  7,  =  330.5  g  cm2,  J2  =  J3  =  5783.0  g  cm2. 

(92a) 

The  initial  conditions  for  the  numerical  models  (IMPACT28  and 
IMPACT35)  are  taken  the  same  as  the  experiment  (see  Sec.  8): 

x0  =  0,  yo  =  0>  Zo  =  0>  u0  =  0,  v0  =  - 1.55  m  s'1, 


Wq  =  —  2.52  ms1, 

»Aio  =  0,  1^20  =  60°,  i//30  =  -95°  ,  to10  =  0, 

co20  =  0.49  s_1,  tn30  =  0.29  s-1.  (92b) 

Substitution  of  the  model  parameters  (92 a)  and  the  initial  condi¬ 
tions  (92 b)  into  IMPACT28  and  IMPACT35  leads  to  the  predic¬ 
tions  of  the  cylinder’s  translation  and  orientation  that  are  com- 


Experlment  3-D  Model 


Fig.  14  Movement  of  cylinder  2  (Z. = 1 2.1 0  cm,p=1.67  g  cm-3)  with  *=-1.00  cm  and  drop  angle 
30  deg  obtained  from  (a)  experiment,  ( b )  3D  IMPACT35  model,  and  (c)  2D  Impact28  model 
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pared  with  the  data  collected  during  Exp-1  at  time  steps  (Fig.  13). 
The  new  3D  model  (IMPACT35)  simulated  trajectory  agrees  well 
with  the  observed  trajectory.  Both  show  the  same  slant-straight 
pattern  and  the  same  travel  time  (1.23  s)  for  the  cylinder  passing 
through  the  water  column.  However,  the  existing  2D  model  (IM- 
PACT28)  has  less  capability  to  predict  the  cylinder’s  movement  in 
the  water  column.  The  travel  time  predicted  by  IMPACT28  is  1.5 
s,  much  more  than  the  observed  value. 

Negative  \  (Nose-Up):  Cylinder  2  (L=  12.10  cm, p 

=  1.67  g  cm-3)  with  ^=-1.00  cm  is  injected  to  the  water  with  the 
drop  angle  30  deg.  The  physical  parameters  of  this  cylinder  are 
given  by 

m  =  254.2  g,  /,=  271.3  gem2,  J2  =  J3  =  3312.6  g  cm2. 

(93a) 

The  initial  conditions  for  the  numerical  models  (IMPACT28  and 
IMPACT35)  are  taken  the  same  as  the  experiment  (see  Sec.  8): 

*o  =  0,  To  =  0.  =  u0  =  0,  v0  =  -  0.75  m  s  , 


w0  =  -  0.67  ms, 

<Aio  =  0,  ifi20  =  24°,  if/30  =  -  96  0  ,  &>i0  =  0, 

&>20  =  -5.08  s"\  &>30  =  0.15  s_1.  (936) 

The  predicted  cylinder’s  translation  and  orientation  are  compared 
with  the  data  collected  during  Exp-1  at  time  steps  (Fig.  14).  The 
new  3D  model  (IMPACT35)  simulated  trajectory  agrees  well  with 
the  observed  trajectory.  Both  show  the  same  flip-spiral  pattern  and 
the  same  travel  time  (1.73  s)  for  the  cylinder  passing  through  the 
water  column.  The  flip  occurs  at  0.11  s  (0.13  s)  after  the  cylinder 
enters  the  water  in  the  experiment  (IMPACT35).  After  the  flip,  the 
cylinder  spirals  down  to  the  bottom.  However,  the  existing  2D 
model  (IMPACT28)  does  not  predict  the  flip-spiral  pattern. 

9.2  Comparison  Using  Exp-2  Data.  After  running  the  two 
models  (IMPACT35  and  IMPACT28)  for  each  gravity  core  re¬ 
gime  [pj(z) , 5(z)]  from  the  initial  conditions  (90),  the  burial 
depths  were  compared  with  measured  burial  depth  data  (Fig.  15). 
As  evident,  IMPACT35  improves  the  prediction  capability.  The 
existing  2D  model  (IMPACT25)  overpredicts  actual  burial  depth 
by  an  order  of  magnitude  on  average.  However,  the  3D  model 
(IMPACT35)  predicts  the  burial  depth  reasonably  well  without 
evident  overprediction.  Since  the  gravity  cores  were  taken  for 
approximately  2  to  3  m  from  the  impact  location,  several  cores 
were  taken  for  each  drop.  This  allowed  an  average  to  be  calcu¬ 
lated  in  order  to  yield  more  accurate  data  for  each  drop. 

10  Conclusions 

1 .  A  3D  model  (IMPACT35)  is  developed  to  predict  the  trans¬ 
lation  and  orientation  of  a  falling  rigid  cylinder  through  air, 
water,  and  sediment.  It  contains  three  components:  triple  co¬ 
ordinate  transform,  cylinder  decomposition,  and  hydrody¬ 
namics  of  a  falling  rigid  object  in  a  single  medium  (air, 
water,  or  sediment)  and  in  multiple  media  (air-water  and 
water-sediment  interfaces). 

2.  Triple  coordinate  transform  is  useful  for  modeling  the  move¬ 
ment  of  a  rigid  body  in  air-water-sediment.  The  body  forces 
(including  buoyancy  force)  and  torques  are  represented  in 
the  E-coordinate  system,  the  hydrodynamic  forces  (such  as 
the  drag  and  lift  forces)  and  torques  are  represented  in  the 
F-coordinate,  and  the  cylinder’s  moments  of  gyration  are 
represented  in  the  M-coordinate.  The  momentum  (moment 
of  momentum)  equation  for  predicting  the  cylinder’s  trans¬ 
lation  velocity  (orientation)  is  represented  in  the 
E-coordinate  (M-coordinate)  system.  Transformations 
among  the  three  coordinate  systems  are  used  to  convert  the 
forcing  terms  into  E-coordinate  (M-coordinate)  for  the  mo¬ 
mentum  (moment  of  momentum)  equation. 


Fig.  15  Comparison  among  observed  and  predicted  burial 
depths 


3.  During  the  penetration,  the  part  that  contacts  the  fluid  (air  or 
water)  is  treated  as  an  equivalent  cylinder  with  the  same 
mass  and  PCOV  location.  The  buoyancy  and  hydrodynamic 
forces  and  torques  are  computed  in  the  equivalent  cylinder. 
The  procedure  developed  for  calculating  external  forcing 
(buoyancy  and  hydrodynamic  forces  and  torques)  for  a 
single  cylinder  is  used  for  the  equivalent  cylinder. 

4.  Impact  force  and  torque  below  the  water-sediment  interface 
are  calculated  on  the  basis  of  the  fact  that  at  the  instance  of 
penetration,  the  sediment  exerts  an  impact  force  only  on  the 
portion  of  the  cylinder’s  surface,  which  moves  towards  the 
sediment.  The  normal  and  tangential  components  of  the  im¬ 
pact  force  are  calculated  separately.  The  normal  component 
is  calculated  using  the  sediment  density  and  shear  strength 
profiles.  The  tangential  component  is  computed  using  the 
friction  law  between  two  solid  bodies  (i.e.,  proportional  to 
the  normal  component).  The  torque  is  easily  obtained  after 
the  impact  force  is  determined. 

5.  The  dynamic  system  for  predicting  trajectory  and  orientation 
of  a  rigid  cylinder  in  air,  water,  and  sediment  are  highly 
nonlinear.  For  example,  the  apparent  torque  in  the  moment 
of  momentum  equation  (20)  (represented  in  the 
F-coordinate)  is  nonlinear.  The  drag  and  lift  forces  are  non¬ 
linear  terms  which  depend  on  the  square  of  the  fluid-to-body 
velocity.  Two  major  assumptions  are  used  to  simplify  the 
system.  First,  the  apparent  torque  is  neglected.  Second,  for 
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the  given  time  step  tn,  the  nonlinear  drag  and  lift  forces  and 
torques  are  linearized  at  any  time  instance  with  temporally 
varying  coefficients  (also  dependent  on  the  fluid-to-cylinder 
velocity)  evaluated  at  the  previous  time  step,  tn_\ .  With  the 
given  cylinder’s  parameters  (translation,  velocity, 
orientation,  and  angular  velocity)  at  the  time 
step  t„,  [x(tn)  ,y(tn)  ,z(tn) ;  u(ta) ,  v(tn),w(tn) ;  i/q  (fn) ,  i/f2(tn) , 
i/r3(rn),&)1(r„),&>2(t„)’<u3(f«)]-  the  model  has  analytical  solu¬ 
tions  at  the  time  step  tn+1.  The  recursive  procedure  is  estab¬ 
lished  to  predict  the  cylinder’s  translation,  velocity,  orienta¬ 
tion,  and  angular  velocity  through  air,  water,  and  sediment 
from  the  initial  conditions.  The  strength  of  such  treatments 
guarantees  the  convergence  of  the  model  integration. 

Since  neglect  of  the  apparent  torque  is  feasible  only  for 
slow  rotation  around  the  cylinder’s  main  axis  (i.e.,  small 
self-spin  angular  velocity  oq),  and  since  local  linearizations 
of  drag  and  lift  forces  and  torques  are  feasible  for  relatively 
small  fluid-to-cylinder  velocity,  the  model  might  not  be  valid 
if  co i  or  the  fluid-to-cylinder  velocity  is  large  such  as  for  fast 
water  entry  and  fast  self  spinning.  A  fully  numerical  calcu¬ 
lation  (rather  than  the  recursive  procedure)  should  be  devel¬ 
oped  for  the  prediction. 

6.  Two  cylinder  drop  experiments  were  conducted  to  evaluate 
the  3D  model.  Model-data  comparison  shows  that  IM- 
PACT35  improves  the  prediction  capability  drastically  with 
an  order  of  error  reduction  in  the  cylinder  burial  depth,  more 
accurate  cylinder  track  (depth  and  orientation)  prediction, 
and  more  accurate  travel  time  of  the  cylinder  through 
air-water-sediment. 
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